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Abstract The adjustment of parameterized gravity waves associated with model convection and finer 
vertical resolution has made possible the generation of the quasi-biennial oscillation (QBO) in two 
Goddard Institute for Space Studies (GISS) models, GISS Middle Atmosphere Global Climate Model III and a 
climate/middle atmosphere version of Model E2. Both extend from the surface to 0.002 hPa, with 2° x 2.5° 
resolution and 1 02 layers. Many realistic features of the QBO are simulated, including magnitude and 
variability of its period and amplitude. The period itself is affected by the magnitude of parameterized 
convective gravity wave momentum fluxes and interactive ozone (which also affects the QBO amplitude 
and variability), among other forcings. Although varying sea surface temperatures affect the parameterized 
momentum fluxes, neither aspect is responsible for the modeled variation in QBO period. Both the 
parameterized and resolved waves act to produce the respective easterly and westerly wind descent, although 
their effect is offset in altitude at each level. The modeled and observed QBO influences on tracers in the 
stratosphere, such as ozone, methane, and water vapor are also discussed. Due to the link between the gravity 
wave parameterization and the models' convection, and the dependence on the ozone field, the models 
may also be used to investigate how the QBO may vary with climate change. 


1. Introduction 

The quasi-biennial oscillation (QBO) is the prime mode of variability in the tropical stratosphere and a large 
contributor to interannual variability throughout the middle atmosphere. It has thus been a long-term goal to 
have models simulate it. Two of the primary forcing functions, mixed Rossby-gravity waves and Kelvin 
waves, are sufficiently large scale that they can be resolved by current climate models if they have sufficient 
vertical resolution. Mixed Rossby-gravity waves can have vertical wavelengths of 3-5 km, and if four points 
are necessary to resolve a wave, resolution on the order of 1 km would be required. This is often achieved in the 
troposphere of global climate models but not necessarily in their stratosphere. But these waves do not seem to 
be sufficient [e.g., Takahashi and Boville, 1992]. 

The additional, smaller-scale waves that seem to be required [e.g., Dunkerton, 1997] are more problematic. 
The very fine vertical and horizontal resolution required for their generation, propagation, and dissipation, the 
hydrostatic dynamics generally employed in climate models, and the coarse parameterizations for potential 
generating mechanisms, such as convection, provide a challenge for global scale models to produce these 
waves and have them propagate and interact properly. 

Hence, with the horizontal grid scales global climate models (GCMs) often use, they have trouble generating 
a realistic QBO without any parameterized aid for the small-scale waves. One exception occurs when the 
inherent subscale dissipation is reduced by up to an order of magnitude [ Takahashi , 1996, 1999], which 
implies that numerical noise may act as a primary (hence unphysical) source or that higher diffusion obstructs 
the normal wave damping mechanisms. In addition, it is recognized that the variability in precipitation- 
induced latent heating is responsible for driving resolved waves in models [e.g., Horinouchi, 2002]. Some 
components of a QBO have been generated with sufficiently fine horizontal and vertical resolution (e.g., 1° 
horizontal resolution, 800 m vertical resolution; [ Hamilton et al., 1999, 2001]) when the model uses a moist 
adiabatic adjustment scheme, but that is thought to overestimate the variability and hence the resolved 
wave generation [Horinouchi 2002; Horinouchi et al., 2003]; indeed, the variability is highly dependent on the 
convection scheme employed [e.g., Ricciardullu and Garcia, 2000]. Even in these cases, much of the 
momentum-flux spectrum likely remains unresolved [Fritts et al., 2006], and a truly direct generation is likely 
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impossible as long as convection is parameterized (otherwise requiring very fine horizontal scales). Recent 
research in this area involves manipulating more sophisticated convection schemes such as Arakawa- 
Schubert so as to increase convective variability for the sake of generating the QBO [ Kawatani et at., 2009]. 
With such changes and as computer power increases, finer resolution models are starting to appear that have 
the capability to generate QBO-like oscillations directly (e.g., 60 km horizontal resolution, 300 m vertical 
resolution; [Watanabe et at., 2008, Kawatani et at., 2011; Evan et at., 201 2]), although many of the problems 
discussed above still remain, resulting in issues such as phase locking to the annual cycle and reduced 
amplitudes in the lower stratosphere [ Kawatani et at., 2011]. 

On the other extreme, modelers can impose a QBO directly by relaxing back to the observed QBO-varying 
stratospheric winds [e.g., Balachandran and Rind, 1995]. This allows one to explore how the QBO winds 
and temperatures can influence other regions. It does not allow one to assess how the QBO could change 
with climate, and there is the possibility that the relaxation process influences the results in ways that differ 
from how the atmosphere would respond were the QBO to be generated more self-consistently. 

An intermediate approach is to use parameterized gravity waves either in a general "nonorographic" 
framework [Scaife et ai., 2000; Giorgetta et at., 2002, 2006] or associated with the expected generating 
mechanism (i.e., convection). The latter approach represents a chief advance in the field, with the 
parameterization of sources becoming process specific. Kim et at. [201 3], using the Met Office Unified Model, 
relate part of their parameterized tropical gravity wave flux to the cloud top momentum flux, as estimated from 
diabatic heating. Richter et at. [2014] with the Community Atmosphere Model relate the gravity wave flux to 
convective heating depth, amplitude, and the mean wind. Schirber et at. [2014], with the ECHAM6 model, relate 
their parameterized momentum flux to convective heating depth, heating rate, and the background wind. 
While the flux magnitude, phase velocity, and wavelength of the parameterized waves are chosen so as to help 
the model simulate at least some aspects of the QBO phenomena, there are some constraints on these choices 
from existing observations [e.g., Alexander et ai., 2010]. 

In contrast to the relaxation approach, the convective gravity wave parameterization effect is on the source 
of the required waves, the propagation and the breaking mechanism, not on the QBO phenomena 
(winds and temperatures) directly. The hope is that this "reduction" in level of control will allow the 
modeled atmosphere a greater freedom to respond interactively, so that the QBO influence is simulated 
more realistically. With the general nonorographic forcing approach, the choice of gravity wave forcing is 
more ad hoc, unconnected to tropospheric sources, often seasonally and annually independent. As 
noted by Kim et ai. [201 3], the generalized nonorographic forcing "results in restricted spatiotemporal 
variations of the parameterized gravity waves (GWs) in models, in spite of a great degree of (real world) 
variability in the GW sources, especially convection in the tropics." Furthermore, there is no indication as to 
how the QBO forcing could change with climate, outside of variations in parameterized propagation 
associated with changes in the intervening winds. These limitations are relaxed when using the forcing 
dependent upon the model's convective activity itself (either parameterized or when generating the waves 
directly [Beres etal., 2005; Richter etal., 2010; Kawatani etai., 2011]). Results do depend on the validity of the 
convective parameterization, whose uncertainties also affect many other phenomena (in addition to the 
QBO and its potential change with climate), e.g., tropical rainfall, Hadley Cell dynamics but for which there 
are many more observable parameters to act as constraints. 

This use of the parameterized gravity wave forcing dependent upon the model convection is employed 
here to generate a QBO in two different GISS models. It is the first of a set of papers that attempt to evaluate 
the effect of climate change on the QBO and, conversely, the effect of the QBO on climate. Results will 
only be as valid as the assumptions noted above, i.e., that the use of parameterized sources does not 
corrupt the atmospheric response elsewhere and that convection (and its change with climate) is handled 
at least somewhat realistically in the model. We will compare with observations when available to assess 
the first aspect. The accuracy of the convective parameterization is in general covered in the source 
references for the models themselves. While this paper is concerned with the generation of the QBO in the 
two models, anticipated subsequent papers will discuss the modeled influence of the QBO on the 
troposphere, stratosphere, and mesosphere in the current climate (Part 2), the effect of climate change on 
the QBO and its influence (Part 3), and the sensitivity of the atmospheric system to joint QBO and solar 
forcing variations (Part 4). 
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Table 1. Vertical Resolution Averaged Over the Layers in Different 
Atmospheric Regions in the Tropics (4°N-4°S) for Both GISS Models, 
Model III and Model E (in Parenthesis), Used in This Paper 3 


Region 

Pressure Levels (hPa) 

Average Thickness (m) 

Troposphere 

Lower 

1 000-800 

142 (140) 

Middle 

800-400 

268 (265) 

Upper 

400-1 00 

357 (353) 

Stratosphere 

Lower 

100-50 

527 (533) 

Middle 

50-30 

707 (688) 


30-10 

1013 (1020) 

Upper 

10-5 

1174 (1171) 


5-1 

1349 (1348) 

Mesosphere 

Lower 

1-0.1 

2901 (2819) 

Middle 

0.1-0.01 

3586 (3537) 

Upper 

0.01-0.002 

4492 (4508) 


a The thicknesses (m) are averaged over the last 50 years of the runs. 


2. Model Description 

The two models used for this study had a 
common origin but have since diverged 
substantially. The GISS Global Climate Middle 
Atmosphere Model III was specifically 
designed to incorporate middle atmosphere 
processes in a climate model, e.g., with finer 
vertical resolution in the stratosphere, a higher 
model top, and a fuller suite of gravity wave 
drag parameterizations. GISS Model E was 
targeted primarily for climate change studies 
with a particular emphasis on tropospheric 
radiative forcing and ocean-atmosphere 
interaction. A discussion of their differences at 
the time was provided by Rind et at. [2007]; the 
subsequent evolution of Model E is provided 
by Schmidt et at. [2014]. A summary of the 
current relationship is provided here. 


Focusing first on the similarities which encompass the very basic structure: the two models have generally 
similar atmospheric radiation schemes based on the correlated K distribution approach for thermal 
radiation gaseous absorption [Lads and Oinas, 1991] and the Single Gauss Point adaptation of the doubling 
and adding equations for shortwave radiation [ Lads and Hansen, 1 974]. Their dynamical numerical 
schemes are quite similar, except for the placement of the gravity wave drag (GWD), as noted below. Both 
models use a mass flux cumulus parameterization originally described by Del Cenio and Yao [1993], 
based on a cloud base neutral buoyancy flux closure with two entraining plumes sharing the mass flux. 
Stratiform clouds are based on a Sundqvist-type prognostic cloud water approach with diagnostic cloud 
fraction [ Del Genio et at., 1996], The boundary layer and surface schemes have as a basic framework, 
the work of Hartke and Rind [1997] and Cheng et at. [2002], while the land surface scheme is derived from 
Friend and Kiang [2005]. 


Flowever, in each subroutine there are now numerous differences. The atmospheric radiation code contains 
much more sophisticated aerosol physics in Model E, with altered absorption coefficients for some gases. The 
details of the convection scheme, i.e., the calculation of updrafts, downdrafts, entrainment rates, convective 
condensate, and physics time step all differ with the result that convective mass fluxes are quite different, 
even when the two models are run at the same resolution. This will become apparent below, when the 
convective gravity wave momentum fluxes for the two models are shown. Similarly, the details of cloud 
microphysics (supersaturation, condensation, and evaporation) and precipitation initiation have diverged, as 
have the boundary layer turbulence formulations and land surface characteristics (soil biophysics and 
vegetative components). For more information concerning these differences, see the discussions in Rind etal. 
[2007] and Schmidt et al. [2014]. 


2.1. GISS Global Climate Middle Atmosphere Model III 

The first model used for these experiments is the GISS Global Climate Middle Atmosphere Model III, run at 
2°x2.5° resolution with 102 layers extending from the surface to 0.002 hPa [Rind et al., 2007]. The vertical 
layering in the tropics for this model (as well as for Model E, whose thicknesses differ due to a slightly 
different temperature structure) is shown in Table 1, where it can be seen that in the troposphere (below 
100 hPa), the resolution is less than 500 m at all levels. In the lower stratosphere (up to 50 hPa), it is about 
500 m. It gradually increases to 1 km by 1 0 hPa, near 1 .5 km at the stratopause. The version used also includes 
the linearized ozone scheme (LINOZ) of McLinden etal. [2000] plus a simple parameterization for tropospheric 
ozone as well as a methane tracer [Rind et al., 2007]. 

The model has parameterized gravity wave momentum fluxes associated with flow over topography, 
convective activity, wind shear, and wind deformation [Rind et al., 1 988, 2007]. The waves "propagate" vertically 
starting from above the source, i.e., in the case of convection, above the top of the convective plume. 
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The propagation is dependent upon background wind conditions by being continually compared against 
ambient saturation levels, with momentum deposition occurring when these are exceeded or as they 
approach critical levels. The wave is assumed to remain saturated and hence its amplitude is assumed to 
remain constant through the "breaking" region (instead of growing exponentially with the density 
decrease, [Lindzen, 1981]), and once above can then resume propagating normally to higher levels to 
deposit additional momentum. All the momentum is deposited at the model top, if not beforehand 
[Lindzen, 1985]; hence, it complies with that aspect of the momentum conservation constraint [Shepherd 
and Shaw, 2004]. (There is a form of Rayleigh friction near the model top, to help counter its influence, 
but its location at 0.002 hPa should minimize its impact on the stratosphere, Shawetal. [2009].) Momentum 
is also conserved in the process as whatever momentum the parameterized waves add to higher 
altitudes is removed from the depth of the source region (so for penetrating convective waves between the 
level of convective origin and the level of convective termination). The kinetic energy lost in the process 
is conserved by being put back as heat locally. To allow the dynamics to fully adjust to the drag, the 
parameterization is incorporated into the numerical solution of the momentum equation (in Model III), 
calculated 30 times per hour. In model E, discussed below, it is imposed as an independent momentum 
forcing calculated in a separate subroutine. 

For the parameterized convective waves, of predominant interest here, the momentum flux is a function of 
the convective mass flux, Brunt-Vaisala frequency at the top of the convective region, wind velocity 
averaged over the convective layers, and horizontal wave number oriented in the direction of the wind 
[Rind et at., 1988, equation (7)]. The horizontal wave number is integrated from 1° resolution to the size of 
the grid box, with a weighting function inversely proportional to the wave number (assuming shorter 
wavelengths influence a more restricted portion of the grid). In practice the wavelengths vary between 
80 km at the highest latitudes and 270 km at low latitudes, the expected order of tropical gravity waves that 
can propagate into the stratosphere. For more details, see Rind et at. [1 988, 2007]. The mass flux in the 
model is strongly related to the depth of penetration, and thus this parameterization is somewhat similar to 
that of the other models discussed in the introduction that use convective sources. 

The actual QBO forcing may be associated with radiative damping of waves as they approach a critical level. We 
incorporated a simple radiative damping parameterization for the parameterized waves dependent upon 
the wave's vertical group velocity in each layer [Rind et at., 1988] and found that it contributed only about -10% 
to the wave forcing (without significantly changing the total forcing or its altitude); the resultant QBO was 
not significantly different. Apparently the model does not care about the exact mechanism of momentum 
deposition. Krismerand Giorgetta [2014] evaluated the damping of resolved waves in the MPI-ESM model and 
showed that radiative damping is only important for large-scale waves. For wave numbers >10, as is the case 
here, horizontal diffusion is the main mechanism of wave damping. For simplicity, the runs were performed 
without this radiative damping effect. 

The standard version of the model did not feature a QBO, and when it was increased in resolution to rx1°, it 
still did not. We did, however, find that the standard model could generate a QBO if the GWD 
parameterization was altered in the following ways: 

1. Phase velocity: The standard model used parameterized convective gravity waves with phase velocities 
±1 0 m s _1 in the direction of the local background wind for shallow convection and additional 
±20 m s _1 and ±40 m s _1 in the direction of the background wind for convection that penetrates above 
400 hPa. For example, if the background wind over the convective altitude is toward the west at —7 m s _1 
for nonpenetrating convection, the parameterized waves would be given phase velocities of —17 m s _1 
(toward the west) and ±3 m s _1 (toward the east). If the convection was penetrating, additional phase 
velocities of -27 m s _1 and ±13 m s _1 would be used (from the ±20 m s _1 component) and —47 m s _1 and 
±33 m s _1 (from the ±40 m s _1 component). Note that these are initial source phase velocities; whether 
the waves propagate into the stratosphere depend on the tropospheric winds. 

The highest phase velocities were instituted to allow the waves to propagate into the mesosphere [e.g., 
Holton, 1983]. Tests indicated that the ±40 m s _1 parameterized waves were not necessary to maintain a 
reasonable mesospheric simulation. Furthermore, their inclusion provided forcing in the upper stratosphere 
that shifted the QBO upward (peak amplitude now near 4hPa) and diminished its amplitude in the middle 
and lower stratosphere. As this was not desirable, they were removed. 
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However, four additional waves (±10 m s _1 and ±20 m s -1 ) orthogonal to the background zonal wind 
were added, providing essentially isotropic source directions, although with varying source phase 
velocities. For the example given above, this would add parameterized waves with meridional components 
of ±10 m s- 1 and —10 m s _1 for nonpenetrating waves and additionally ±20 m s _1 and —20 m s _1 for 
penetrating waves. (Note that opposite azimuth waves do not "cancel," because each is "launched," and 
their propagation and breaking level likely differ). This was done to better capture the filtering effect of 
upper tropospheric winds; as discussed below, this filtering leads to a seasonality of the effective 
stratospheric forcing and appears at least partly responsible for many aspects of the QBO, including its 
variability. We note that the use of essentially isotropic source directions is consistent with the approach of 
Choi and Chun [2011] from mesoscale numerical simulations, while the incorporation of the background 
wind in the determination of the phase velocity is also generally consistent with their use of the moving 
speed of the convective source as determined by the background wind (one difference: Choi and Chun 
[2011] use the wind at 700 hPa, while here it is averaged over the convective altitude range). 

Models attempting to simulate the QBO often use a greater spectrum of phase velocities in excess of the 
eight employed here. Tests were done with additional values (e.g., increments of ±2.5 m s _1 , ±5 m s _1 , 
etc.). They had no obvious effect on the QBO simulation, partly because many of them were absorbed 
within the troposphere and also because the eight that are being used cover a varying range already, 
given the temporal and spatial variations in the background wind magnitude upon which they are 
superimposed. In addition, we also experimented with up to 18 different azimuths but again found little 
improvement in the simulated QBO. Ultimately we opted for the minimum number of both phase 
velocities and azimuths deemed sufficient, keeping in mind the prospective use of this parameterization 
in long climate change simulations. Note that the momentum flux deposition ultimately calculated is 
applied to both the zonal and meridional wind field. 

2. Convective gravity wave flux magnitude: The "constant" coefficient in Rind et al. [1988, equation (7)], 
relating the convective gravity wave momentum flux to the model's convective mass flux, has two 
components: a scaling term and a factor accounting dimensionally for area weighting and time 
integration of the convection. The scaling term was increased compared with that employed in the 1 02 
layer model of Rind et al. [2007] such that it produces momentum fluxes some 3 times larger in the 
convective regions, amounting to 2 times larger on the global average. These values are for the experiment 
with an embedded scaling coefficient of 0.143 (henceforth, "MU scale"). Conceptually the MU scale can be 
thought of as implying either the fraction of the grid box undergoing convection and/or the intermittency 
with which the convective waves are being generated during the 1 h physics time step. In the experiments 
described here, that coefficient was altered within the range of 0.13 to 0.15. We note that a similar 
"efficiency" coefficient relating to the fractional area of convection was employed by Richter et al. [201 4], 
which they increased from 0.1 to 0.55 to provide sufficient forcing for their QBO simulation. 

These changes have the effect of increasing the number of parameterized convective gravity wave groups 
from six to eight, with greater azimuthal coverage and increasing the amplitude of their momentum flux. 

In no other respect was the earlier 2x2.5, 102 layer model changed, e.g., the other parameterized gravity 
waves were not altered. It appears as if the increased azimuthal coverage was necessary in our model to allow 
for both propagation of the waves to the required altitudes and for the "critical level" deposition mechanism 
to be sufficiently activated. The increase in momentum flux was required for wave breaking to occur at 
appropriate levels. Given that the background wind velocity is often from the southeast or northeast, the 
orthogonal component is adding an additional zonal component as well. The larger amplitudes are also 
consistent with what seems to be required in models, as discussed below. 

The absolute value of the parameterized moist convective momentum fluxes (the sum of the absolute value 
of the fluxes over all four azimuths) at the tropospheric source altitude are shown for the two solstice seasons 
and on the annual average in Figure 1 (left) (averaged over the ensemble of runs for Model III). Shown 
separately are the absolute values with ±10 m/s and ±20 m/s including values both aligned and orthogonal to 
the background wind. The full zonal average tropical magnitudes for the sum of these contributors (1.5-2 mPa) 
fall in the range of the typical nonorographic forcing used in some models (1-5 mPa reported by Geller et al. 
[201 3]) and are in general agreement with the forcing used by Kim et al. [201 3] associated with their model's 
convection. Peak tropical values in Model III of 4-6 mPa are in agreement with that shown at 1 00 hPa in July 
from gravity-wave-resolving models [Soto et al., 2009; Alexander et al., 201 0]. Hence, compared with other 
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MC Momentum Fluxes 


Model III 


Model E 


10 m/s 
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Max 3.7 0.26 Global 


DJF 



Max 1.8 0.17 Global 



Max 4.3 0.29 Global 


JJA 



Max 5.4 0.28 Global 


Max 2.3 0.19 Global 



60 

30 
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-30 
-60 
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.0 .4 .8 
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1.6 




2.0 


2.4 2.8 


Figure 1 . Absolute value of the convective gravity wave momentum fluxes averaged over the ensemble of runs (Model III) or the length of the run (Model E) for the 
two solstice seasons and the annual average. Results are shown for the parameterized fluxes associated with phase velocities ±10 m/s and ±20 m/s along and 
orthogonal to the background wind. The global average and maximum value are shown for each figure. Shown are results for (left) Model III and (right) Model E. 


models, the previous values used in this model were too small in amplitude. The geographic and seasonal 
variations are associated with the model's varying moist convection. As a matter of note, this altered convective 
parameterization (and others) were also included in the 2x2.5, 53 layer model (with the same model top, 
hence decreased vertical resolution); while a QBO was obtained, it did not descend sufficiently into the lower 
stratosphere compared with observations and thus would be expected to have reduced influence in the 
troposphere [Garfmkel and Hartmann, 2011]. The increased (by a factor of 2) vertical resolution of -500 m 
available with the 1 02 layer model (Table 1 ) is in accordance with what has often been deemed necessary in other 
QBO model studies [e.g., Takahashi, 1 996; Horinouchi and Yoden, 1 998; Hamilton et al., 2001 ; Richter et at., 2014]. 

2.2. GISS Model E2 

In the course of this work, for comparison purposes, we also created an equivalent version of the GISS Model E2 
([ Schmidt et al., 2014], henceforth simply Model E), starting from the version used in the Coupled Model 
Intercomparison Project Phase 5 (CMIP5) database, but with 102 vertical layers, the top raised to 0.002 hPa (Table 1, 
thickness values in parenthesis) and the complete suite of parameterized gravity wave drag formulations included. 
Several of the parameterized gravity wave momentum fluxes, for flow over topography and model-resolved 
deformation, were already in Model E, and their wave parameters were retained; the additional gravity wave 
formulations for shear and convection utilized the formulations from Model III, although the convective MU scale 
was calibrated to produce the approximate observed QBO period. The formulation for the drag near the top of the 
model in the CMIP5 version was retained when the top was moved to the higher levels; the consequences of this 
will be discussed below. In all other respects, the CMIP5 version formulation and parameterizations were utilized. 

As detailed above, for many of their physical parameterizations, Model III and Model E are quite different. 
Nevertheless, the convective gravity wave drag parameterizations used for Model III also produced a QBO in 
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Figure 2. Tropical zonal winds as a function of altitude and time averaged between 4°N and 4°S for the simulation of (top) Model III (with MU scale = 0.1 43), (middle) 
Model E, and (bottom) between 5°N and 5°S in observations (from ERA40, Uppata eta!., [2005]). Results are shown for the common time frame available from the ERA40. 


Model E (it did not have one previously), with a minimal amount of tuning of the MU scale. In this case, one 
simulation was performed, from 1880 to 2010. It included observed sea surface temperatures and sea ice, 
changes in solar irradiance, and volcanic aerosols, but unlike Model III, also included changes in greenhouse 
gases and aerosols. The run was made with a calibrated value of the MU scale of 0.1 8 providing close to the 
observed QBO period. (Note that Model E has different convective mass fluxes from Model III, so a comparison of 
MU scales is not overly relevant.) 

The absolute value of the convective gravity wave momentum flux calculated in this model is shown in Figure 1 
(right); the time-averaged tropical values in Model Eof ~1 mPa are similar to those obtained from satellites at 20 km 
(as summarized by Geller et at. [201 3]). As shown in the figure, the peak tropical values of the moist convective 
momentum fluxes are somewhat lower than in Model III, with the global mean convective momentum fluxes 
about 20% lower. This version of Model E had decadal variations in atmospheric composition using precalculated 
stratospheric and tropospheric ozone fields so as to provide the radiative forcing for the transient CMIP5 
simulations (as opposed to calculating the ozone fields and their perturbations via LINOZ as in Model III), and the 
prescribed values did not contain an ozone QBO. (Other versions of Model E do include fully interactive ozone 
chemistry [Shindell et at., 201 3].) The impacts and consequences of this aspect will be discussed below, in section 4. 

3. Simulation of the QBO 

3.1. Mean Characteristics 

The GISS Middle Atmosphere Global Climate Model III with the modifications indicated above was integrated 
with the observed sea surface temperatures, solar forcing, and volcanic aerosols (hence "natural forcings") 
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Monthly Zonal Mean Wind (U) 

Averages 


32 hPa MU Scale = .143 
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Figure 3. Tropical zonal winds as a function of time at 62 hPa and 32 hPa from (top) Model III (MU scale = 0.1 43) and at 
70 hPa and 30 hPa from (bottom) observations. Observations are from the tropical wind analyses at the Tropical Islands 
(FU-Berlin data set) produced by combining the observations of the radiosonde stations at Canton Island (3°S, 1 72°W), Gan/ 
Maledive Islands (1°S, 73°E), and Singapore (1°N, 104°E) (Free University of Berlin, n.d, http://www.geo.fu-berlin.de/en/met/ 
ag/strat/products/qbo/). 


from 1951 to 2004. The simulations were performed six different times, for a total of 324 years of results. 
Each simulation differed only in the precise choice of MU scale, varying from 0.13 to 0.15. Shown in 
Figure 2 (top) are the zonal winds averaged from 4°N to 4°S as a function of altitude and time, with 
the MU scale = 0.1 43 in Model III; the period of the wind oscillation and its peak power vary with the 
choice of this scale, as detailed below. Also shown are the results for Model E (Figure 2, middle) and 
observations from the ERA40 data [Uppala et at., 2005] averaged from 5 N to 5°S, featuring a similar 
time frame as the model (Figure 2, bottom). In both models the QBO dominates the variability from 
about 1 00 hPa to 5 hPa, with the semiannual oscillation taking over in the upper stratosphere, consistent 
with observations. 

Pawson et at. [1993] and Giorgetta et at. [2006] noted a somewhat faster and more regular downward 
propagation of the westerly wind QBO phase (hence, less time spent at any level in the easterly phase, 
at least at some levels), with the westerly wind phase propagating without loss from approximately 
1 0 hPa to 50 hPa. The models' relative phase propagation can be explored by assessing the length of time 
spent in the easterly wind phase versus the westerly wind phase at a given pressure level for each cycle. 
The results are given below, utilizing for the observations the FU-Berlin data set of tropical island 
observations (Free University of Berlin, n.d, http://www.geo.fu-berlin.de/en/met/ag/strat/products/qbo/, 
henceforth, "Tropical Islands"). 
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Figure 4. Tropical zonal winds as a function of time at 62 hPa and 32 hPa in Model E. 


First, at ~30 hPa (mean/standard deviation in months): 

Tropical Islands (Westerly 13.4/2.7; Easterly 14.3/2.8); 

Model III, MU scale = 0.1 43 (Westerly 16.0/2.8; Easterly 1 1.6/1 .0); 

Model E for the same 55 years, (Westerly 14.7/2.0; Easterly 13.6/2.1). 

Then at ~60 hPa: 

Tropical Islands (Westerly 17.0/3.9; Easterly 10.4/3.6); 

Model III (Westerly 19.4/2.3; Easterly 9.1/1 .4); 

Model E (Westerly 10.0/3.0; Easterly 17.4/4.5). 

The values indicate that at 30 hPa there are no significant differences between the length of the easterly and 
westerly phases in either the models or the observations. At 60 hPa, Model E spends relatively too little time 
in the westerly wind phase. 

Pawson et al. [1 993] and Ciorgetta et al. [2006] also noted a somewhat stronger intensity of the easterly winds 
that diminishes somewhat toward lower levels. The greater intensity of the easterly wind relative to the 
westerly wind can be seen more clearly for models and observations in Figures 3 and 4. Shown in Figure 3 
(top) are the results for Model III with MU scale = 0.143 at 62 hPa and 32 hPa. The winds oscillate between 
±1 5 m/s at 62 hPa and from +20 to —30 m/s at 32 hPa, in approximate agreement with observations shown in 
Figure 3 (bottom) (whose 70 hPa lower level is at a slightly greater pressure level). An equivalent picture 
for Model E is provided in Figure 4, in which it can be seen that the QBO wind amplitudes are somewhat 
smaller than for Model III (with that particular MU scale). 

In order to assess all of these aspects in a clearer framework, we have produced a "composite" QBO picture 
for the two models and the ERA40 observations. This was done by using the amplitude of the zonal winds at 
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Figure 5. Composite QBO for (top) Model III, MU Scale = 0.1 43, (middle) 
Model E, and (bottom) ERA40 data. Results are for the years 1951-2004 
for the models and 1958-2001 for the observations, averaged between 4°N 
and 4°S. See text for details. 


30 hPa to define the phase for each 
cycle (0° being the transition from 
easterly to westerly, 90° maximum 
westerly wind, etc.). Each QBO cycle was 
sampled 12 times (every 30°). The 
dates for each phase were then used to 
generate an average "snapshot" of 
the QBO tropical zonal wind structure 
(4°N-4°S). The results are presented in 
Figure 5 for the two models and 
observations as a function of the 
phase angle; while this angle does not 
correspond to a fixed month 
throughout the various QBO cycles, the 
0° phase is often around January at 
30 hPa in both the models and 
observations. With an approximate 
28 month period covering the 360°, that 
would then imply that Winter is roughly 
centered around 0° and 150°, Spring 
around 60° and 210°, Summer at 90° 
and 240°, and Fall at 1 20° and 270°. The 
general correspondence between the 
models and ERA40 data is apparent, as 
well as some differences, primarily that 
the models' westerly winds are a bit 
stronger than the observed, as depicted 
in the ERA40 reanalysis, the models' 
easterly winds somewhat less "thick" in 
altitude and their descent a little slower. 
With respect to the Model E wind 
duration discrepancy presented above, 
westerly winds at 60 hPa in Model E 
arrive at approximately the 120° phase, 
whereas in both Model III and ERA40 
they arrive about 30° earlier, indicating 
that their descent from 30 hPa in Model 
E is not sufficiently fast. We return to 
this result and the somewhat lower 
QBO wind amplitudes for Model E 
compared with Model III in section 4, 
along with the influence of ozone. 


Krismer et al. [2013] in their simulation 
with the Max Planck Institute Earth 
System Model noted that at 5 hPa, the 
onset of the westerly wind regime 

clustered in spring and fall due to the coupling of the QBO with the semiannual wind oscillation, which at that 
altitude (—34 km in the tropics) is westerly in those seasons; the primary onset occurred in May. Using the 
approximate season/phase values given above, with respect to Figure 5, this would be at approximately 220° 
(at 5 hPa). As noted in their paper and apparent in the figure, in the ERA40 data, the onset arrives in June 
(-230°). For Model III, it arrives in April (-210°) and for Model E in February-March (-180°). 

Krismer et al. [2013] also found that the 5 hPa westerly winds are stronger in (Northern Fiemisphere) spring 
than in fall in the ERA40 observations and their model and propagate downward faster. The authors attribute 
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Figure 6. Periodograms of the zonal wind at various pressure levels from the MU scale = 0.1 43 run (left) for Model III, and 
(right) for Model E. The peak period is indicated by dashed lines and on the top of each panel. The purple line at 0.0075 hPa 
for Model III is from separate simulations with reduced Rayleigh friction near the model top. 


this to the seasonal modulation of the QBO forcing, i.e., enhanced wave filtering by the lower stratosphere 
westerly winds in fall and winter, compared with spring and summer. This would imply that the westerly wind 
downward slope with time would be greater at 60° and 210° than at 120° and 270°, which is true in both 
the models and ERAA40 observations. The spring westerlies then descend from 5 to 50 hPa in a little over a 
year [ Krismer et at., 201 3], and so both westerly and easterly wind regimes occur primarily during Northern 
Hemisphere spring at ~50 hPa [Dunkerton, 1 990; Pawson et a!., 1 993], In the ERA40 data, the westerly 
wind regime appears to maximize from winter to spring (140°-210°), with the easterly wind regime 
maximizing from winter to early spring (330°-30°). Corresponding model values at 50 hPa are, for both 
Model III and Model E, westerly maximum in winter to early spring (120°-180°), easterly maximum in late 
winter to early spring (0-30°). Given that the westerly wind initiation began a little earlier in the models 
than ERA40 at 5 hPa, the westerly wind descent is approximately the same in the model and observations. 
The easterly wind descent looks a little slower in the models. An exception to all this occurs at times of 
the presence of substantial volcanic aerosols, when the resultant tropical heating in the lower stratosphere 
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Figure 7. As in Figure 6 except for observations, from the (left) Tropical Islands and the (right) ERA-Interim reanalysis. 

can lead to an extended period of westerly winds at that altitude (shown for example in Figure 2 from the 
Pinatubo aerosol in 1993). 

3.2. Period of the QBO 

The periodogram of the zonal wind at these and also higher altitudes (lower pressures) is provided in 
Figure 6, for Model III with MU scale of 0.143 (left) and for Model E (right). For comparison we show the 
periodogram at many of these same altitudes from observations in Figure 7, for both the Tropical Islands 
and the ERA40 Interim reanalysis [ Dee etal., 2011], The numbers in each panel indicate the peak period, and, 
in Figure 7, its associated power. As expected, peak periods vary from close to 28 months at the lower 
altitudes, to the semiannual oscillation in the upper stratosphere (which is a resolved feature of the model 
without the parameterized convective gravity waves). These values are also presented in Table 2 for 30 hPa 
and 60 hPa. 

The mean period can also be calculated based on the phase of the wave, e.g., the length of time spent before 
a particular phase repeats (such as transitioning from westerly to easterly wind). Adding the average easterly 
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and westerly wind durations given in the previous subsection produces the average period by this phase 
transition approach. The results can be compared to the value given in Table 2 from the calculated 
periodogram, and as can be seen, the two approaches give very similar numbers. 

Model E even features a noticeable mesospheric QBO (MQBO) as well as the expected more predominant 
semiannual oscillation at 80 km, and the MQBO increases with altitude to this level in the general manner 
reviewed by Baldwin et al. [2001]. This is in contrast to Model III, which has no substantial MQBO peak. 

A difference between the two models is that the Rayleigh-type friction, used near the model top to reduce 
its influence (the gravity wave drag reaching that level being not sufficient) is twice as large in Model III as 
it is in Model E, so the stronger drag may have acted to inhibit the MQBO appearance. To test this 
hypothesis, a subset of runs was carried out with Model III, utilizing the same (reduced) mesospheric drag 
as in Model E. The MQBO was now somewhat more coherent, as indicated by the purple line in Figure 6. 
In terms of integrated power for the MQBO, the value in Model III was 40% of that in Model E, while in 
the new simulation, it has increased to 80%. Hence, most of the difference is due to this difference in 
mesospheric drag. 

The period of the QBO oscillation varies with the magnitude of the parameterized momentum flux, as with 
stronger momentum deposition, the respective easterly and westerly winds are drawn downward faster [e.g., 
Giorgetta et al., 2006]. This is quantified in Table 2, where the peak period and associated power are compared 
for the runs with different global momentum forcing and observations at 30 hPa and 60 hPa. Note that the 
power at the peak is greater in the ensemble than in any of the individual runs, as the signal comes through 
more clearly. That is also true for the observations; comparison of the Tropical Islands result for the full time 
period (1953-201 1, as in Table 2) with values for the somewhat shorter time period of 1979-201 1 as in Figure 7, 
similar to that of the ERA40 Interim analysis, shows that the power increased 40-60% at the different altitudes 
with the longer record (hence explaining most of the difference between the Tropical Islands and ERA40 
Interim). Thus, for comparison's sake, we show for Model E the results for a similar time period as in Model III and 
the observations; the results indicate a general agreement between the models and observations. The 
difference in QBO period between the Tropical Islands and ERA40 Interim is entirely the result of the different 
years utilized; the ERA40 Interim project employed the tropical island winds in their reanalysis. 

To counter the "duration" effect on peak power, we also show in Table 2 the "integrated" power, integrated 
between the peak period and the periods where it has fallen to 'h the peak power. The power is now much 
less dependent on record length, although the ensemble value is reduced due to its combination of different 
inherent periods. Both Model E and the individual experiments in Model III once again have power values 
similar to the observed, at both 30 and 60 hPa. 

While much of the variation among the different experiments can be related to statistical sampling, the 
integrated power assessment at 30 hPa shows a systematic decrease as the momentum flux increases. 
Stronger momentum fluxes would be expected to break at lower altitudes, and given that the phenomenon 
is initiated via interaction with the semiannual oscillation in the upper stratosphere, the momentum breaking 
may be less effective when it is further removed in altitude. The integrated power assessment at 60 hPa 
shows a less systematic version of the same relationship. 

The momentum flux is not the only variable that affects the mean period. Simulations were performed with 
various climate forcings and boundary conditions. The use of solar cycle variability systematically increased 
the period of the QBO by 3/4 of a month; a full exploration of this effect will be undertaken in paper IV of this 
series. The use of volcanic aerosols also systematically increased the mean period by ~3/4 of a month. 
Volcanic aerosols have the effect of warming the tropical lower stratosphere, thus extending a westerly wind 
component [see, for example, Hamilton, 2002]. The results presented previously include these aspects. Ozone 
influence on the QBO period will be discussed, in conjunction with stratospheric tracers, in section 4. 

3.3. Variability of the QBO Period 

Also shown in Table 2 are measures of the variability of the QBO period within each simulation. Again this is 
calculated in two different ways: the standard deviation via the phase transition from easterly to westerly 
winds and via the width of the half-peak power from the periodograms. Considering first the phase 
transitions (the entry in the standard deviations column), it can be seen that from this method of 
evaluation, Model III does provide the approximate variability. However, Model E variability in period 
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appears to be a little too small at 30 hPa; again, we return to the various aspects of Model E's simulation of 
its QBO period in section 4. 

Another view can be obtained by comparing the width of the half-power around the peak on the periodograms. 
This is included in Table 2 in the last column. The value is affected by record length in the sense that 
the shorter the record length, the greater the value; when the Tropical Island data are rerun for the same 
years as the ERA 40 Interim data, its half-power width is essentially the same (as opposed to the 88% 
decrease shown in Table 2). The greater record length is also why the "ensemble" value is so small. We thus 
again use the same time frame for Model E and Model III, and their values are close to those of theTropical 
Islands at both 30 hPa and 60 hPa. Hence, from this perspective, Model E's variability in period is now 
appropriate; that is true from both approaches for Model III. 

The half-power spread differs considerably from the phase method in the ordering of the period's variability 
among the different experiments. Effectively, it uses all phases of the wave, rather than just the transition 
phase, which may account for the difference, but it may simply be that the differences themselves are not 
significant. Assuming that the different estimates of the variability are members of the same population set, use 
of the t test suggests that at 30 hPa, with the "phase" approach, the variability is significantly greater (95% level) 
with the lowest momentum forcing (the MU scale = 0.1 3 experiment). The simulations with the greatest 
momentum fluxes (MU scales = 0.1485 and 0.15) have the smallest variability with the half-power approach. 
Both results fit conceptually with the concept that some models that exhibit too small variability of the QBO 
period, or values phase locked to the annual or biennial cycle, may have too strong forcing, by resolved waves 
as well as by parameterized ones [e.g., Hamilton et al., 1 999; Takahashi, 1 999; Kawatani et at., 2011]. 

Given the sensitivity of the period to the magnitude of gravity wave momentum flux forcing, it is reasonable to 
suppose that the observed variation in QBO period is directly related at least in part to the variation in 
gravity wave forcing [e.g.. Getter et at., 1997]. The experiments listed in Table 2 can establish a quantitative 
relationship, in this case for the difference in period between the different experiments. Using the tropical 
parameterized momentum fluxes (3°N-3°S, rather than the global values shown in Table 2) originating at the 
top of the convective plume, we find the following; 

1. y= 55 — 16.3 lx; correlation coefficient of —0.56 at 60 hPa 

2. y = 58 — 1 7.79x; correlation coefficient of —0.65 at 30 hPa 

with y the period in months, and x the tropical momentum flux in mPa. 

We can perform a similar linear regression using the varying tropical momentum fluxes and the varying 
periods for each of the QBO cycles during a particular run. The procedure was carried out for Model III, (MU 
scale = 0.143) and for Model E, both at 60 hPa. The results are as follows: 

1. y = 48 — 1 2.1 2x; correlation coefficient of —0.33 for Model III at 60 hPa 

2. y= 42 — 16.41x; correlation coefficient of —0.27 for Model Eat60hPa. 

Hence, between 1 0 and 30% of the variance in the period is related to the variations in tropical momentum 
fluxes in these models from these approaches. The greater momentum forcing draws the easterly and 
westerly regimes downward more quickly. 

To carry this one step farther, a prime initiator of varying convection in the model relates to the variation 
in tropical sea surface temperatures, of which the El Nino-Southern Oscillation (ENSO) phenomena is the 
leading mode of variability. We compared the ENSO 3.4 Index with the Model Ill's parameterized 
convective momentum flux for both the nonpenetrating ( V — 10 m/s) waves and the penetrating convective 
waves ( V - 20m/s), in the MU scale = 0.143 simulation. The correlation with the nonpenetrating 
parameterized forcing was small (r = 0.2 1 ), but it was higher, as expected, with the penetrating convective 
forcing (r=0.48). From this perspective, in this model, 'A of the variance in our higher-phase velocity 
parameterized waves is ENSO related; given the previous quantitative relationships, this would be 
associated with only about 1/20 or about 5% of the variation in QBO period. The effect would likely be too 
small to be seen. Maury et al. [2012] noted that their model's resolved equatorial waves were significantly 
affected by ENSO variations, and Yang et at. [201 2] briefly note that stratospheric tropical waves resolvable 
in the ERA40 Interim data also have an ENSO variation. Relating such variations to the QBO period in 
observations is difficult, because the biggest ENSO events have tended to occur around the times of 
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volcanoes, which can affect the QBO period; and given the duration of ENSO events, 1 to 2 years, the QBO 
can, and has, often spanned a transition from positive to negative SSTs or the reverse. 

In fact, other modeling efforts have produced realistic variations in QBO period while using climatological sea 
surface temperatures and even constant tropospheric gravity wave sources. The seasonally varying winds 
near the tropopause impose a filter on parameterized waves entering the stratosphere, and it is this 
seasonality that is thought to be the strongest producer of variability of QBO forcing [Krismer et al., 2013; 
Schirber et al., 2014]. However, even if true, varying convective GWD sources may have an effect on other 
aspects, such as the seasonal cycle of phase propagation discussed above. 

To quantify the impact of varying gravity wave sources or sea surface temperatures in a GISS model, we 
performed several additional experiments. In the first, the convective gravity wave source momentum, its 
phase velocity, its geographic location, and its level of launching were saved for each time step for 1 year 
(1950) with Model E, and the same forcing was then utilized each year, on the same date. The model was 
allowed to filter the waves differently, due to its normal background wind changes, so the stratospheric 
impact was allowed to change. The variability of the QBO period was then compared to the same 
simulation years with its normal varying parameterized sources. The results, given in Table 2 as "Model E, 
Const Mom Flux," show that neither the variability of period defined by the standard deviation from its 
phase or the half-power spread has been reduced; the differences compared to Model E run normally are 
not significant. 

The parameterized forcing used here is only one component of the total QBO forcing that will likely have 
interannual variations, as the resolved large-scale wave forcing will also likely be altered due to variations in 
sea surface temperatures/convection and probably other effects as well. These were all free to change in the 
experiment and may have contributed to some of the resultant variation in period Therefore, one additional 
experiment utilized the same specified forcing from the previous run but this time with climatological sea 
surface temperatures and sea ice (1950-2000), hence no interannual variations. The results presented in 
Table 2 as "Model E, Const Mom Flux/Clim. SSTs" again show little systematic change in the measures of 
variability. We conclude from these various perspectives that while the magnitude of the momentum flux 
source does affect the QBO period, in practice, other factors dominate the variability in QBO period. There 
was also little change in the seasonal cycle of phase propagation. 

3.4. Forcing of the QBO 

We next explore the comparative influence of the model's resolved and parameterized tropical waves. The 
models have the horizontal and vertical resolution to resolve mixed Rossby gravity waves (wave number 3-5, 
4-6 day period) and Kelvin waves (wave number 1-3, 10-20 day period), and to some extent inertia gravity 
waves (periods <1 day, wavelengths <1000 km), but, as noted above, that was insufficient to generate a QBO 
without the additional smaller scale parameterized convective gravity waves. The first question to address is 
how well the model is generating the "resolvable" waves or whether the parameterized waves have to 
compensate for these as well. 

Shown in Figures 8 and 9 are the space-time spectral analyses for the tropical zonal and meridional winds at 
50 hPa from 1°N to 1°S in Model III (MU scale = 0.143) and for the equivalent years of Model E, calculated using 
the Maximum Entropy Method for short time periods (i.e., each month) [Hayashi and Golder, 1980; Hayashi, 
1982]. The monthly values were then used to create a long-term (54 years) average for each of the models. The 
power in the U wind for eastward moving waves (i.e., positive periods) at the lowest wave numbers and longest 
periods represent Kelvin waves, and the power in the V wind for westward moving waves (i.e., negative periods) 
at more intermediate wave numbers represent mixed Rossby gravity waves, as well as some Rossby wave 
energy [e.g., Yang et a!., 2012]. While there are some differences between the models, overall their simulation of 
these tropical waves is similar, both qualitatively and quantitatively. 

These results can be compared with those of Yang et al. [201 2] from the ERA40 Interim data (averaged over 
15°N to 15°S, compared with 1°N to 1°S in the models). For Kelvin waves, for example, wave 2 at ~15day 
period, the ERA40 data and both models show power in the U wind of about 0.07 m 2 s~ 2 (the plotted data for 
U in that paper are at 20 mb, compared with 50 mb for the models, but Yang et al. [2012] show that there is 
little change in power between those levels). Peak Kelvin wave power is about 0.1 6 m 2 s~ 2 in the ERA40 data, 
slightly less than the value in Model III and slightly more than in Model E. Overall, the power in the U wind at 
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Figure 8. Space-time wave power spectrum for the zonal and meridional winds in the tropics (1°N-rS) at 50 hPa in Model 
III averaged over 54 years. Results for Model III are from the simulation with MU scale = 0.143. Minimum and maximum 
values are shown in the upper right hand corner. To calculate the energy in (m/s) 2 for an individual wave number/period, 
divide the value in the figure by 1 0 times the period (for the U wind), and 1 00 times the period (for the V wind). 


50 mb, for eastward moving waves with wave numbers between 2 and 1 0 and periods between 2 and 30 days 
is somewhat less than 2 m 2 s~ 2 in the ERA40 data, while it is 1 .5 m 2 s~ 2 in Model III and 1 .4 m 2 s~ 2 in Model E. 

For mixed Rossby gravity waves, for example, wave number 4 at 5 day period, the ERA40 data indicate a power 
in the V wind at 30 mb of 0.03 m 2 s~ 2 ; values for similar waves in Model III and Model E at 50 mb are about 0.1 2 
and 0.09 m 2 s~ 2 , respectively (again the ERA40 data show little change between these altitudes). Peak power 
in the mixed Rossby gravity wave range is about 0.15 m 2 s~ 2 in the ERA40 data, slightly higher than that in 
both models. Overall, the power in the V wind at 50 mb, for westward moving waves with wave numbers 
between 2 and 1 0 and periods between 2 and 30 days is about 1 m 2 s 2 in the ERA40 data, 1 .6 m 2 s~ 2 in Model III, 
and 1.5 m 2 s~ 2 in Model E. The results suggest that the models are simulating these equatorial trapped tropical 
waves with reasonable power; whatever differences that exist with the ERA40 data may be at least partly 
associated with the different latitudinal domains utilized for the calculation. At least in this respect, the models 
may therefore be expected to force the mean flow in a manner not unlike that which occurs in the real world. 

As noted in section 1, it is now recognized that resolved waves are not sufficient in most models to generate a 
QBO or in the assessment from observations, either [e.g., Gray and Pule, 1989; Dunkerton, 1997; Horinouchi, 
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2002], although the amount of additional smaller-wavelength wave forcing appears to be model dependent 
[Horinouchi et al., 2003; Scinocca and McFarlane, 2004]. In particular, both GISS models use penetrative 
convection schemes, including downdrafts that typically do not overestimate tropical precipitation variability 
and, hence, by current model resolution standards, would be expected to require more rather than less 
parameterized waves (although as noted by Scinocca and McFarlane [2004], the large-scale precipitation 
scheme also influences tropical precipitation variability). 

To present the average resolved and parameterized forcing effect on the zonal wind, we make use of the 
approach employed to generate the "composite" picture of the zonal wind in Figure 5. Using the same dates 
for the respective phases, we computed the average zonal wind forcing from the parameterized waves 
(Figure 10, top row), from the resolved waves (via the Eliassen-Palm (EP) flux divergence, Figure 10, second 
row), associated (transformed) advection (Figure 10, third row), and the zonal wind change occurring at that 
time (Figure 1 0, bottom row) for both models. The parameterized moist convective wave drag effect is similar 
in both models, and, as expected, it produces forcing where the appropriate phase velocity waves are 
trapped by the equivalent zonal winds (compare Figure 1 0 with Figure 5). It thus strengthens the zonal wind 
at the altitude at which it is already strong. Flowever, there is also a strong degree of forcing slightly below 
that level, which imparts a tendency for descent. The easterly wind acceleration is greater than the westerly 


RIND ETAL. 


©2014. American Geophysical Union. All Rights Reserved. 


18 




^AGU Journal of Geophysical Research: Atmospheres 


1 0.1 002/2014JD021 678 


4°N— 4°s Avg Composite du/dt Changes 

Model-Ill MUs=. 143 


-17.8 12.4 


MC Drags 


ModelE (last 55 yrs) 



.001 


60 120 180 240 300 360 

-5.7 39.7 



EP Fix Diverg 

r\r\ i 


-6.6 27.5 



0 60 120 180 240 300 360 


0 60 120 180 240 300 360 


-15 


-10 


(10 6 m/ s/s) 


0 


10 


15 


Figure 10 . Composite contributions of their respective effects on the zonal wind (4°N-4°S) of the (top row) parameterized con- 
vective gravity waves, (second row) resolved waves, (third row) advective effect, and its (bottom row) total change for (left) 
Model III with MU Scale = 0.143 and (right) Model E. These correspond to the composite zonal wind results shown in Figure 5. 


wind acceleration in both models, presumably because the parameterized phase velocity, influenced by the 
background easterly winds in the tropics, has a predominantly westward component and thus westward 
momentum flux. 

The resolved waves are in general providing forcing opposite to the sign of the local zonal wind, hence 
instigating its change. They therefore are contributing substantially to the descent of the zonal winds. There 
is an inclination for the EP flux divergence and transformed advection terms to offset one another to some 
extent, as expected from nonacceleration theorem considerations, a result that can also be seen in the 
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Figure 11 . Time varying ozone at 62 hPa and 32 hPaforthe (top) MU scale = 0.1 43 run of Model III and the (bottom) periodogram. 


corresponding breakdown shown by Ciorgetta et al. [2006] from the MAECHAM5 model. Comparing the 
respective forcings to the sign of the total change that is actually occurring (bottom row) shows that at each 
level the parameterized waves are contributing strongly to the upper portion of the total change (i.e., the zonal 
wind descent), while the resolved waves are forcing the middle and lower portions, with a somewhat more 
modest magnitude. Overall, therefore, they are both contributing, though in different ways. Giorgetta et al. 
[2006] with the MAECHAM5 mode found the parameterized and resolved waves to be contributing with 
somewhat similar magnitude. 

Shown in both Figures 5 and 10 are the mesospheric QBO and the contributions to it for Model E. The 
"yellow" bar near the model top at 0.003 hPa in Figure 5 (middle panel) represents the westerly wind phase 
of the mesospheric QBO, centered around 180 phase; this is approximately a 180° phase difference relative 
to the maximum easterly winds at 40 hPa (near 0° phase). Such a phase difference between these two 
levels has been reported from H RDI data [Burrage et al., 1996; Baldwin et al., 2001]. As shown in Figure 10, 
both the parameterized and resolved waves are forcing maximum westerly winds in the mesosphere near 
180 phase; while the resolved wave effect is stronger, it is offset to a good extent by the transformed 
advection. These results agree with the theoretical explanation for the mesospheric QBO as resulting from 
a combination of selective filtering of propagating small-scale gravity waves and wave breaking at upper 
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Figure 12. Tropical (top) ozone anomaly, (middle) methane anomaly, and (bottom) water vapor anomaly as a function of 
altitude and time for a portion of the simulation, with MU scale = 0.1 43 in Model III. 


levels [Baldwin et at., 2001]. The forcings are somewhat similar for Model III, although as noted above, this 
version did not feature a mesospheric QBO, a result that is due mostly to the stronger Rayleigh friction used 
near the model top. 

4. The QBO in Tropical Stratospheric Tracers 

As Model III calculates its ozone field, we can see the QBO impact on this variable as well. Shown in Figure 1 1 
is the time varying field at 62 hPa and 32 hPa, along with the spectrum for the MU scale = 0.1 43, both 
depicting the QBO. The tropical ozone anomaly as a function of altitude for MU scale = 0.1 43 is presented in 
Figure 12 (top) for a portion of the run. It is very similar to that derived from Stratospheric Aerosol and Gas 
Experiment data as shown by Randel and Wu [1 996] as well as Tian et al. [2006], including the abrupt phase 
change at ~28km, an effect which was also simulated in the CCM coupled climate-chemistry model [Tian 
et al., 2006], As discussed by Baldwin et al. [2001] this level separates the region of the atmosphere below 
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Figure 13 . Composite QBO cycle for modeled (top left) ozone, (top right) methane, (bottom left) water vapor, and (bottom 
right) temperature. 


where ozone can be treated as a long-lived tracer from the region above, where its short lifetime makes it 
responsive to chemical and hence thermal perturbations. The model's results arise from both the QBO- 
induced temperature change impact on ozone photochemistry and transport change. Note that this model 
does not include QBO variations in NO y transport, one proposed mechanism for the ozone QBO; Butchart 
et al. [2003] also found that it was not needed in their model to produce this effect, while Tian et al. [2006] 
found that it was operative primarily around 35 km, hence above the level of phase reversal. 

Tropical easterly or westerly wind fields require from thermal wind constraints opposite changes in the 
latitudinal temperature gradient at levels below, accomplished by vertical motions. As the QBO zonal wind 
changes with altitude (e.g., Figure 2), the circulation must also reverse at that altitude with the vertical 
velocity acting on the vertical ozone gradient. The reversal of the circulation with altitude during Northern 
Hemisphere winter is consistent with the residual circulation change shown by Randel et al. [1 999], and the 
implied meridional wind changes are generally consistent with those derived from observations by Hitchman 
and Huesmann [2009]. 

The residual circulation change also affects long-lived species transported within, and to a lesser extent into, 
the stratosphere. A methane tracer was turned on during these experiments in Model III and is converted to 
water vapor in the upper stratosphere (for more information on it, see Rind et al. [2007]). The results for a 
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portion of the simulation are shown in Figure 12 (middle). The variation is a transport effect, and the modeled 
QBO effect is weak or nonexistent at the lowest stratospheric levels, as methane has only a small vertical 
gradient in the tropics from the tropopause to ~30 hPa altitude [e.g., Baldwin et a!., 2001]. 

The QBO would also be expected to alter the "tape recorder" effect of ascending water vapor anomalies in 
the tropical lower stratosphere due to the seasonal cycle of tropopause temperature. In Model III, the vertical 
velocity in this region averages about 03 mm/s, similar to those derived from the Flalogen Occultation 
Experiment water vapor data [Schoeberl etal., 2008]. The QBO-induced vertical velocity anomaly in the model 
peaks at about Vi that magnitude. The specific humidity anomaly arising in response is shown in Figure 12 
(bottom) for a portion of the record, and the QBO signal induces a perturbation peak of roughly 25% 
background values. Also visible are stronger water vapor intrusions associated with volcanic warming of the 
tropopause due to El Chichon and Mt. Pinatubo. 

Finally, we present in Figure 1 3 the QBO composite variation of these fields in the stratosphere, along with 
the composite temperature variation, using the same dates established from the zonal wind variation at 
30 hPa for Figure 5. At levels below 1 0 hPa, the ozone anomaly is in phase with the specific humidity 
and temperature anomalies, and out of phase with that for methane; both ozone and specific humidity 
increase with altitude in the tropical stratosphere, while methane decreases with altitude, so these phase 
relationships emphasize the importance of vertical transports for all three species. Where the momentum flux is 
causing an easterly QBO wind, say at 0° phase angle in the lower stratosphere, it requires from the thermal wind 
relationship a decreased latitudinal temperature gradient relative to the subtropics at levels below, achieved by 
forcing an upward velocity in the tropical lower stratosphere. This results in colder tropical temperatures 
relative to the subtropics, while the upward transports produce reduced values of ozone and water vapor, with 
increased methane. The inverse applies when the momentum transport produces westerly QBO winds. 

At levels above 1 0 hPa, ozone now becomes more influenced by chemical/thermal anomalies associated with 
the QBO. For example, when momentum transports produce a QBO westerly wind field, the increased 
latitudinal temperature gradient necessitated by thermal wind constraints is achieved by downward velocity 
in the tropics. This results in warming. The warmer temperatures inhibit net ozone production, resulting in 
reduced ozone. Given their vertical gradients, it also produces increased water vapor and decreased 
methane. Ozone anomalies are thus now in phase with methane anomalies and out of phase with those for 
water vapor and temperature. In both regions, the ozone changes modify the temperature changes via 
absorption of shortwave and longwave radiation, amplifying it when in phase at levels below 1 0 hPa, and 
diminishing it when out of phase, above. 

As noted above, the version of Model E used here employed an input ozone field (as opposed to calculating it 
interactively), and it did not contain a signature of the QBO in the ozone fields. To test the impact of this 
difference, we reran Model III with that same ozone input field (and hence without the model-produced ozone 
QBO); the results, shown in Table 2 as "0.143 w/o LINOZ" indicated that the period of the QBO was about 
2.5 months shorter without the coupling (compare with the 0.143 run), similar to the differences found by 
Butchart et al. [2003] and Tian et al. [2006] in coupled climate-chemistry models. The variability of the QBO 
period was also reduced by -25% (Table 2), especially with respect to the calculation based on phase at 30 hPa; 
this might then be the reason for the lower value in Model E. We note also that were the background ozone 
profiles to change with climate, so would the ozone QBO, and hence the QBO amplitude and variability. 

In this study, and also that of Tian et al. [2006], the changes in ozone affected both the absorption of solar 
radiation and of outgoing infrared radiation thus affecting the modeled temperature QBO. Flad this diabatic 
effect of chemistry coupling been included in Model E, then to produce a "proper" period of -28 months would 
have required an increase in convective momentum flux of -20% (Table 2, assuming a similar response of the 
two models), bringing the momentum flux for Model E into better agreement with that for Model III. 

Tian etal. [2006] also showed that the ozone coupling increased the amplitude of the QBO, and that was seen 
in these experiments as well, with the zonal wind peaks about 25% larger at 30 and 60 hPa when ozone was 
responding; it is also apparent in the difference in the power numbers in Table 2 for this experiment. 
Previously, using a 2D model, Dingmin et al. [1996] had calculated an ozone-induced 25% effect on the QBO 
temperature response and 10% on the zonal wind. The ozone influence, combined with the reduced 
parameterized convective momentum fluxes, are likely responsible for the somewhat reduced QBO zonal 
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wind amplitudes seen in Model E. They do not, however, by themselves explain the differences in length of 
easterly and westerly wind phases seen between this model and observations at 60 hPa; the relative easterly 
and westerly wind durations at this level were not affected in the "rerun" of Model III without the ozone QBO. 
However, Model E's QBO westerly winds were weak at 60 hPa and were the convective gravity wave 
momentum fluxes to be recalibrated in conjunction with the coupled ozone impact, they would be 
presumably be stronger and perhaps longer lasting. 

5. Discussion 

A simulation of the QBO has long been a goal of modeling groups, and the results show that a QBO-like 
oscillation has been achieved in two different GISS models. The use of the model's convection to produce 
parameterized gravity waves that induce a QBO is the key feature of such experiments, differing from the 
"nonorographic" forcing in use by some modeling groups. To the degree it is successful, it does not "prove" 
that the parameterization is duplicating the actual sources that help force the observed effect. It only 
demonstrates that it can be done this way. There is no guarantee that the smaller-scale waves deemed 
necessary from the observations actually arise from convection, although we know that convection plays a 
significant role in the generation of tropical waves, and the high-frequency waves utilized here are known to 
be associated with deep convection [e.g., Alexander and Holton, 1997]. 

The exact importance of convection to the generation of the QBO is still a matter of investigation. Diabatic 
heating associated with convection can be seen to directly drive slower-phase velocity equatorial waves in 
the low-to-middle troposphere, while in the upper troposphere and lower stratosphere, equatorial waves 
with higher phase velocities propagate apart from this heating, as so-called dry or free modes [e.g., Kiladis 
et a!., 2009]. Also, see Yao and Jablonowski [201 3] for a simulation without convection capable of producing 
the relevant tropical waves. 

By allowing the parameterized waves to propagate vertically during each hour, we are assuming they have 
fast vertical group velocities [e.g., Dunkerton, 1 997], We use certain phase velocities and orientations of the 
parameterized waves, along with a particular scheme for parameterized wave-mean flow interaction, while 
the real world mechanisms are undoubtedly much more complicated. However, the phase velocities and 
wavelengths used are consistent with those derived from observations [e.g., Bergman and Salby, 1994; Sato 
et al., 1 995]. The ability of this simplified approach is to some extent the result of the various calibration 
experiments performed to achieve it, but it may also indicate that the system is robust enough to produce 
the observed response if given the chance (e.g., sufficient vertical resolution and sufficient upward 
momentum fluxes, here related to convection). At least in these GISS models, the details of the respective 
(differing) model physics does not appear to be a determining factor. 

6. Conclusion 

Following are the primary results from this study: 

1 . The QBO has been generated in two GISS models, at 2° x 2.5°, 1 02 layers extending from the surface up 
to 0.002 hPa. 

2. It was done primarily by increasing the number of parameterized convective gravity waves from 6 to 8, 
with greater azimuthal coverage, increasing the amplitude of their momentum flux by 3 times in 
convective regions and increasing the vertical resolution to -500 m in the lower stratosphere. 

3. The modeled QBOs share with the observed numerous characteristics, such as the proper period and 
variability, the onset of the westerly wind regime clustering in spring and (secondarily) fall at 5 hPa, 
faster and more regular downward propagation of the westerly wind phase, and a somewhat stronger 
intensity of the easterly winds. 

4. The primary deficiency in both models is that the westerly winds are stronger than observed in the midstra- 
tosphere, and the easterly winds are somewhat less deep than observed, and their descent a little slower. 

5. The QBO period varies with the magnitude of the parameterized convective momentum flux, as with stron- 
ger momentum flux deposition, the respective easterly and westerly winds are drawn down faster. 

6. Experiments with constant parameterized convective momentum flux and constant sea surface tem- 
peratures indicate that, in these models, variations in these parameters do not increase the variability 
of the period of the modeled QBO. 
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7. The power in the models' resolved Kelvin and mixed Rossby gravity waves compares well with 
observed values. 

8. Within each level of descent, the parameterized waves are contributing strongly to the upper portion of 
the total change, while the resolved waves are forcing the middle and lower portions, with a somewhat 
more modest magnitude. 

9. A mesospheric QBO is also being forced by both the parameterized and resolved waves. 

10. QBO variations in ozone, methane, and water vapor arise and agree with observed variations. At levels 
below 1 0 hPa, where transports and vertical gradients dominate the changes, the ozone anomaly is in 
phase with the specific humidity and temperature anomalies, and out of phase with methane. 

11. Above lOhPa, where ozone is more responsive to chemical/thermal anomalies, ozone anomalies are 
now in phase with methane anomalies, and out of phase with those for water vapor and temperature. 

12. Use of interactive ozone, and thus the generation of an ozone QBO, lengthens the QBO period by 
2.5 months, and increases both its zonal wind amplitude and the variability of its period by about 
25%. Changes in the background ozone profile with climate may thus affect the QBO. 

While the use of parameterized convective gravity waves is undoubtedly simplistic compared with real world 
forcing, the values chosen for their parameters are broadly consistent with observations, and the similarities of 
the simulated and the real world QBO suggest that the approach should be useful in investigating how the QBO 
influences other regions, how it responds to climate forcings, and how climate change will affect the QBO 
influence on the stratosphere and troposphere. Some previous studies have concluded that an intensification 
of the stratospheric Brewer-Dobson circulation is likely as climate warms [e.g., Rind et at., 1 990, 1 998; Butchart 
et at., 2006; Garda and Randet, 2008]. The stronger upwelling would then produce a tendency for a slower 
downward progression of the wind structures, resulting in an increase in period [e.g., Kawatani et at., 2011], or 
even a substantial weakening in its downward penetration [Kawatani and Hamiton, 2013; Rind et at., 2013]. 
However, with the potential for an increase in penetrative convection as climate warms due to the greater static 
energy available, there would be an opposing tendency of greater convective gravity wave momentum flux 
and a decrease in period. How these two effects compare obviously depends on their respective strengths 
[Watanabe and Kawatani, 201 2]. These are the subjects of subsequent papers in this series. 
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